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In  the  Compressed  Sensing  (CS)  Imaging  with  Wide  FOV  and  Dynamic  Magnification 
Project  Report  we  present  the  development  of  diverse  imaging  systems  based  on  Compressive 
Sensing  or  Compressive  Sampling  (CS)  method.  In  particular,  we  demonstrate  three  different 
implementations  of  the  CS-based  imaging  systems,  including:  (a)  a  single  pixel  camera,  (b) 
a  CS-based  multispectral  imaging  system,  and  (c)  a  CS-based  optical  sectioning  microscope 
(CSM). 

1  Single  Pixel  Camera 

Most,  optical  images  are  sparse  in  some  transformed  domains,  such  in  the  Fourier  or  Wavelet 
domain.  To  exploit  the  sparse  information  in  the  transformed  domain,  we  used  CS  measure¬ 
ment  patterns  to  modulate  the  intensity  of  optical  images.  Based  on  the  CS  measurement 
results,  the  original  optical  image  can  be  reconstructed  by  solving  a  minimization  problem. 
To  experimentally  realize  the  CS  measurement  process,  we  utilized  a  Digital  Micromirror 
Device  (DMD)  to  implement  the  CS  measurement  patterns.  The  core  component  of  the 
DMD  is  a  7G8(V)?1024(H)  aluminum  micromirror  array.  Each  mirror  can  be  electrically 
driven  to  be  turned  on  or  off.  The  DMD  is  installed  in  the  image  plane  of  an  objective  lens. 
When  CS  measurement  patterns  are  being  displayed  by  the  DMD,  in  its  reflection  direction, 
we  can  see  the  intensity  of  optical  images  modulated  by  the  active  DMD  pattern.  The  CS 
measurements  are  collected  into  a  single  pixel  detector  by  a  focusing  lens.  Since  only  a  single 
pixel  detector  is  used  to  acquire  the  complete  information  of  the  optical  image,  this  CS-based 
imaging  system  design  can  also  be  considered  as  a  single  pixel  camera.  Figure  1  shows  the 
schematic  drawing  of  the  single  pixel  camera. 


Figure  1:  Schematic  drawing  of  the  DMD-based  single  pixel  camera. 

Figure  2(a)  shows  the  experimental  setup  of  the  single  pixel  camera.  Figure  2(b),  (c), 
and  (d)  show  one  of  the  CS  measurement  patterns  used  in  our  experimental  setup,  its  DMD 
implementation,  and  an  image-modulation  result  by  the  DMD  pattern,  respectively. 
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Figure  2:  (a)  The  experimental  realization  of  a  single  pixel  camera,  (b)  A  CS  measurement 
pattern,  (c)  The  DMD  implementation  of  the  CS  measurement  pattern,  (d)  An  image- 
modulation  result  using  the  DMD  implemented  CS  measurement  pattern. 


Figure  3  shows  an  image  reconstruction  result  obtained  from  our  experimental  setup. 
Fig  3(a)  shows  a  CCD  image  of  the  imaging  target,  which  is  a  piece  of  paper  printed  with 
two  characters.  Figures  3(b)  shows  the  reconstructed  image  of  the  imaging  target.  In  this 
experiment,  we  used  a  variable  density  sampling  method  in  the  Iladamard  domain  to  realize 
the  CS  measurement  process.  The  reconstructed  image  has  a  pixel  size  of  2567256,  and  it 
was  reconstructed  using  6554  CS  measurement  results.  In  this  result,  we  can  see  that  the 
important  information  of  the  imaging  target  was  successfully  reconstructed. 


¥) 


Figure  3:  (a)  A  CCD  image  of  the  imaging  target,  (b)  An  image  reconstruction  result 
generated  using  our  experimental  setup. 

The  single  pixel  camera  can  also  be  used  to  realize  passive  CS  imaging  in  the  near  infrared 
(NIIl)  spectrum.  Figure  4  shows  the  reconstructed  optical  images  in  the  wavelengths  of  1200 
nm  (Fig.  4(b))  and  1450  nin  (Fig.  4(c)).  Figure  4(a)  shows  the  imaging  target,  which  is 
a  USAF-1951  resolution  target.  In  this  experiment,  NIR  Light  Emitting  Diodes  (LEDs)  of 
emission  wavelengths  of  1200  nm  and  1450  nm  were  used  for  illumination  and  a  single  pixel 
Germanium  (Ge)  detector  was  used  to  collect  the  CS  measurements. 
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(a)  (b)  (c) 

Figure  4:  (a)  Imaging  target,  (b)  Optical  image  reconstructed  in  the  wavelength  of  1200 
nm.  (c)  Optical  image  reconstructed  in  the  wavelength  of  1450  run. 

2  Compressive  Sampling  Multispectral  Imaging  Sys¬ 
tem 


The  hardware  architecture  of  the  single  pixel  camera  can  be  conveniently  extended  to  build 
a  CS-based  Multi-Speetral  Imaging  (CS-MSI)  system.  Conventional  MSI  systems  usually 
use  sensor-arrays  or  single  pixel  detectors  plus  mechanical  scanning  instruments  to  acquire 
optical  images  at  multiple  spectral  channels.  In  a  CS-MSI  system,  expensive  sensor-arrays 
or  complicated  mechanical  scanning  instruments  can  be  avoided  and  only  single  pixel  dev 
tcctors  arc  needed  to  acquire  images  at  multiple  spectral  channels.  This  means  the  eco¬ 
nomic/ technical  investments  on  MSI  systems  can  be  significantly  reduced. 


Figure  5:  Schematic  drawing  of  a  4-channel  CS-MSI  system. 

Figure  5  shows  the  schematic  drawing  of  the  CS-based  MSI  (CS-MSI)  system.  In  this 
design,  the  DMD-reflected  image  irradiance  was  split  into  four  spectral  channels  by  beam¬ 
splitters.  Different  spectral  filters  were  installed  at  the  end  of  the  spectral  channels.  Three 
visible  wavelengt.li  bandpass  filters,  centered  at  650  nm,  550  nm  and  450  nm  were  used  at 
channels-1,  2,  and  3,  respectively.  At  channel-4,  a  1200-nm  long-pass  filter  was  used  to  iso 
late  the  NIR  spectral  information.  The  transmission  bandwidth  of  the  450-nm  filter  is  40 
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nm  whereas  the  550-nm  and  650- nm  spectral  filters  had  a  transmission  bandwidth  of  10  run. 
In  the  NIR  spectral  channel  (Channel-4),  a  single  pixel  detector  based  on  Ge  material  was 
used.  In  the  visible  spectral  channels,  photomultiplier  tubes  (PMT)  were  used.  The  imaging 
target  was  a  color  checker  card  (X-Rite  Mini  ColorChecker).  A  150- W  quartz  halogen  lamp 
(Dolan-Jenner  DC-950H  DC-Regulated  Fiber  Optic  Illuminator)  was  used  to  illuminate  the 
imaging  target,  along  with  a  150-W  Xenon  lamp  (Newport  Research  Are  Lamp  Source) 
to  enhance  the  illumination  in  the  UV  and  NIR  spectrums  w'hen  needed.  Figure  6  shows 
an  image  reconstruction  results  obtained  using  our  setup  (the  top  row  images),  which  are 
compared  with  CCD  images  captured  using  the  same  spectral  filters  as  the  ones  installed  in 
the  CS-MSI  system  (the  bottom  row  images).  CCD  images  are  scaled  and  resized  to  match 
the  reconstructed  images.  No  brightness/contrast  adjustment  was  applied  to  the  displayed 
images. 


Figure  6:  Comparison  between  the  spectral  images  reconstructed  from  the  CS-MSI  testbed 
(top  row)  and  CCD  spectral  images  captured  with  the  same  filters  (bottom  row),  (a)  (b), 
(c)  and  (d)  are  spectral  images  reconstructed  from  the  650-nin,  550-nm,  450-mn  and  NIR 
spectral  channels  respectively,  (f),  (g),  (h)  and  (i)  are  CCD  images  captured  with  the  650- 
nm,  550-nm,  450-nm  and  the  NIR  spectral  filters  respectively,  (e)  Color  image  reconstruction 
result,  which  is  generated  by  filling  the  visible  spectral  images  into  an  RGB  data  format,  (j) 
RGB  synthesis  of  the  CCD  visible  spectral  images. 

Figures  6  (a),  (b),  (c)  and  (d)  are  spectral  images  reconstructed  in  the  650-nrn,  550-nm, 
450-nm  and  NIR  spectral  channels.  Figures  6  (f),  (g),  (h)  and  (i)  arc  CCD  images  captured 
with  the  650-nm,  550-nm,  450-nm  and  the  1200-nm  long-pass  filters.  In  Fig.  6  (i),  we  do 
not  see  any  useful  information  because  the  CCD  camera  used  in  this  experiment  was  not 
sensitive  to  wavelengths  in  the  NIR  spectrum.  To  capture  these  CCD  spectral  images,  the 
DMD  in  the  CS-MSI  system  was  replaced  with  an  exposed  CCD  sensor  array  (AVT  PIKE 
F-100C).  Spectral  filters  were  placed  directly  in  front  of  the  exposed  CCD.  Figure  6  (c) 
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is  generated  by  filling  the  images  reconstructed  from  the  visible  spectral  channels  into  an 
RGB  data  format.  Figure  6  (j)  is  the  corresponding  RGB  synthesis  of  CCD  spectral  images. 
The  reconstructed  images  shown  in  Fig.  6  have  a  pixel  size  of  2567256  (65536  pixels).  The 
Hadamard-space  variable-density  sampling  method  was  used  and  the  down-sampling  ratio 
was  10%  (6534  CS  measurements). 

3  Compressive  Sampling  Optical  Sectioning  Microscope 

In  this  effort,  we  also  implemented  the  CS  method  in  building  optical  sectioning  for  con- 
focal  microscopy  imaging  systems.  In  conventional  optical  sectioning  microscopes,  optical 
sections  are  captured  using  a  single  pixel  detector,  such  as  a  Photo  Multiplier  Tube  (PMT), 
in  a  point-by-point  fashion.  To  acquire  the  complete  information  from  the  imaging  scene, 
usually  a  2-axis  scanning  process  is  needed.  This  scanning  process  can  be  time  consuming. 
To  reduce  the  image  acquisition  time,  a  novel  design  called  Programmable  Array  Micro¬ 
scope  (PAM)  was  proposed.  PAM  systems  use  a  I)MD  to  implement  multi-pinhole  masks 
in  conjugate  image  planes  of  a  confocal  microscope.  Compared  to  conventional  optical  sec¬ 
tioning  microscopes,  the  image  acquisition  speed  of  PAM  systems  is  much  faster,  because 
in  this  case,  the  information  from  multiple  points  in  the  imaging  scene  is  acquired  at  the 
same  time.  In  PAM  systems,  CCD  cameras  are  used  to  acquire  the  confocal  information 
from  the  imaging  target.  In  low  light  or  high  speed  imaging  settings,  CCD  cameras  pose 
high  economic  and  technical  challenges  to  the  PAM  system.  If  CS  measurement  patterns 
are  implemented  with  the  PAM  system,  confocal  images  can  be  captured  by  single  pixel  de¬ 
tectors  and  thus  the  technical/economic  investments  on  CCD  cameras  in  PAM  systems  can 
be  reduced.  PAM  systems  implemented  with  CS  measurement  patterns  can  be  considered 
as  Compressive  Sensing  Microscopes  or  CSMs.  Figure  7  shows  the  schematic  drawing  of  the 
CSM  system. 

Figure  8  shows  our  experimental  realization  of  the  CSM  system  and  one  of  the  image 
reconstruction  results  obtained  using  this  experimental  setup.  The  imaging  target  used  in 
this  experiment  is  a  slide  of  skin  section.  The  objective  lens  used  in  this  setup  is  a  60x 
Olympus  lens,  with  a  Numeric  Aperture  (NA)  of  0.85.  A  532  nm  solid  state  laser  (30  mW) 
was  used  for  illumination.  We  implemented  the  Hadamard  space  variable  density  sampling 
method  in  order  to  realize  the  CS  measurement  process  and  we  also  used  6534  measurement 
patterns  in  this  experiment.  The  reconstructed  image  had  a  pixel  size  of  1287128.  The 
excitation  filter  used  in  this  experiment  had  a  center  wavelength  of  543  nm.  Its  transmis¬ 
sion  bandwidth  was  30  nm.  The  emission  filter  had  a  center  wavelength  of  610  nm  and  a 
transmission  bandwidth  of  75  nm.  In  the  experiment,  we  used  Modified  Scrambled-Block 
Hadamard  Ensembles  (MSBHEs)  as  measurement  patterns.  The  measurement  patterns  had 
a  pixel  size  of  1287128.  The  block  size  (BS)  of  the  MSBHE  patterns  was  16. 

Figure  9  shows  a  wide  field  multispect.ral  microscopic  image  reconstruction  result  (2 
spectral  channels)  obtained  using  the  experimental  setup.  Figure  9  (a)  is  the  imaging  target 
captured  by  a  CCD  camera.  The  imaging  target  was  a  15-?m  fluorescent  beads  slide.  Figure  9 
(b)  shows  the  multispectral  image  reconstruction  result  captured  by  our  experimental  setup. 
Figure  9  (c)  shows  the  green  spectral  image  of  the  imaging  target,  which  was  captured  under 
the  illumination  of  a  532  nm  solid  state  laser.  The  fluorescent  filter  set  used  to  capture  this 
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Figure  7:  Schematic  drawing  of  the  CSM  system. 


spectral  image  was:  excitation  filter  (center  wavelength  (mil) /bandwidth  (mil)):  543/30. 
emission  filter:  010/75.  Figure  9  (d)  shows  the  red  spectral  image,  which  was  captured 
with  a  632  nm  when  a  He-Ne  laser  was  used.  The  fluorescent  filter  set  used  to  capture  this 
spectral  image  was:  excitation  filter  (center  wavelength  (nm)/handwidth  (nm)):  020/00, 
emission  filter:  700/75.  The  filter  sets  are  purchased  from  Chroma. 

We  also  studied  the  optical  sectioning  imaging  performance  of  the  CSM  system.  Figure 
10  shows  6  optical  sections  of  a  reflective  imaging  target  (a  microcircuit  chip)  reconstructed 
using  our  experimental  setup.  In  this  experiment,  we  used  a  white  light  source  (Dolan-Jenner 
DC-950H  DC-Regulated  Fiber  Optic  Illuminator)  for  illumination. 

Figure  10(a)  shows  the  imaging  target.  Figures  10(6)  —  (^7)  show  the  reconstructed  optical 
sections.  The  adjacent  optical  sections  are  1  gw  away  from  each  other.  We  can  six;  the 
reconstructed  information  of  the  imaging  target  changes  as  the  optical  section  moves  in  the 
depth  direction.  Figure  11  show's  another  set  of  optical  sectioning  imaging  results  captured 
using  the  experimental  setup.  In  this  case,  a  biological  specimen  (pollen  grain)  was  used  as 
the  imaging  target.  A  532  nm  solid  state  laser  was  used  for  illumination.  The  excitation  filter 
used  in  this  experiment  had  a  center  wavelength  of  543  11m,  and  a  bandwidth  30nm.  The 
emission  filter  had  a  center  wavelength  of  GlOnra,  and  a  bandwidth  of  75/?m.  The  image  in 
the  top-left  corner  is  the  imaging  target,  which  is  a  wide-held  image  captured  using  a  CCD 
camera  and  white  light  illumination.  The  right  side  images  are  the  reconstructed  optical 
sections  of  the  imaging  target.  The  adjacent  optical  sections  are  1/zm  away  from  each  other 
in  the  depth  direction.  From  this  result,  we  can  see  that  as  the  optical  section  moves  in  the 
depth  direction,  different  information  of  the  pollen  grain  specimen  in  the  depth  direction 
was  reconstructed. 
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Figure  8:  (a)  Experimental  setup  of  the  CSM  system,  (b)  Biological  sample  (skin  section) 
reconstruction  result,  from  a  compressive  sampling  PAM  setup. 

4  Reconstruction  Algorithm 

By  exploiting  the  sparsity  of  natural  images,  compressed  sensing  has  shown  that  it  is  feasible 
to  acquire  and  reconstruct  natural  images  from  a  limited  number  of  linear  projection  mea¬ 
surements  at  sub-Nyquist  sampling  rates  [1,  2).  A  key  to  the  success  of  CS  is  the  design  of 
the  measurement  ensemble,  which  is  based  on  the  evaluation  of  the  incoherence  between  the 
measurement  ensemble  and  the  sparsity  basis.  The  incoherence  property  requires  that  tin* 
components  of  the  measurement  ensemble  should  have  dense  representations  in  the  sparsity 
basis.  Due  to  the  large  scale  nature  of  images,  the  generation  of  t  he  measurement  ensemble 
should  be  both  computationally  efficient  and  memory  efficient.  Furthermore,  the  sampling 
scheme  should  enable  fast  algorithms  to  perform  image  reconstruction. 

Measurement  matrices  where  each  entry  is  an  independent  and  identically  distributed 
(i.  i.  d.)  random  variable  obeying  Gaussian  or  Bernoulli  distribution  have  been  proposed 
for  compressed  image  sampling  [1,  2],  Recently,  it  has  been  shown  that  the  performance  of 
CS  sampling  can  be  improved  if  the  random  measurement  matrices  are  suitably  optimized 

[3] .  These  methods  lead  to  optimized  but  unstructured  measurement  matrices  and  thus  large 
memory  space  and  computation  demanding  resources  are  needed,  making  them  prohibitively 
expensive  for  implementation.  A  more  desirable  way  to  obtain  linear  measurements  is  by 
incoherent,  sampling  in  a  transform  domain  that  is  equipped  with  fast  transform  algorithms 

[4] .  Measurement  ensembles  in  the  transform  domain  that  enable  fast  computations  include 
partial  Fourier  ensemble,  scrambled  Fourier  ensemble  (SFE),  scrambled  block  Hadamard 
ensemble  (SBHE)  and  Noiselets.  These  have  been  shown  to  obtain  good  performance  in 
CS  applications  [2,  4,  5].  With  the  exception  of  signal  sparsity  on  a  given  basis,  however, 
the  formulation  of  these  sampling  approaches  does  not  exploit  any  a  priori  information  of 
natural  images  that  could  lead  to  improved  CS  performance  [6,  7]. 
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(c)  (d) 


Figure  9:  (a)  Imaging  target  used  in  this  experiment,  which  is  a  fluorescent  bead  specimen. 
The  diameter  of  the  beads  is  15-?m.  (b)  Multispectral  image-reconstruction  result  obtained 
from  the  experimental  setup  (2  spectral  channels),  (c)  Spectral  image  reconstructed  when 
a  532  nm  solid  state  laser  wras  used  to  offer  illumination.  The  fluorescent  filter  set  used  to 
capture  this  spectral  image  is:  excitation  filter  (center  w’avelength  (nm) /bandwidth  (nm)): 
543/30,  emission  filter:  610/75.  (d)  Spectral  image  reconstructed  when  a  632  nm  He-Ne 
laser  was  used  for  illumination.  The  fluorescent  filter  set  used  to  capture  t  his  spectral  image 
is  (excitation  filter:  620/60,  emission  filter:  700/75). 

We  developed  a  method  to  effectively  exploit  such  a  prion  information  in  the  design 
of  efficient  CS  measurement  ensemble  by  using  the  inherent  statistical  distributions  that 
natural  images  exhibit  in  the  sparse  wavelet  domain.  A  novel  family  of  variable  density 
sampling  patterns  are  designed  in  the  frequency  transform  domain,  which  is  applicable  to  the 
Discrete  Fourier  Transform  (DFT),  Discrete  Cosine  Transform  (DCT)  and  ordered  discrete 
Hadamard  transform  (DHT)  [8].  To  design  the  proposed  measurement  matrices,  the  widely 
used  generalized  Gaussian  distribution  (GGD)  model  is  adopted  to  statistically  describe 
the  distribution  of  image  w'avelet  coefficients.  The  method  of  incoherent  Fourier  sampling 
of  subband  wavelets  proposed  in  [4]  provides  an  efficient  w'ay  to  acquire  sparse  subband 
wavelets  in  the  DFT  domain.  The  Fourier  coefficients  in  the  band  where  significant  energy 
of  the  wavelet  exists  are  sampled  randomly  to  minimize  the  coherence  between  the  sparse 
wavelets  and  the  measurement  Fourier  atoms.  Based  on  such  incoherent  Fourier  subband 
sampling  strategy,  w'e  derive  a  variable  density  sampling  function  p{rn,  ri)  in  t  he  frequency 
transform  domain  according  to  the  adopted  statistical  wavelet  model.  Here  p{m,  n)  indicates 
the  probability  that  the  (m,n)th  coefficient  is  sampled.  The  design  of  the  variable  density 
sampling  function  is  then  further  improved  so  that  not  only  the  distribution  of  significant 
wavelet  coefficients  are  considered,  but  also  their  rapid  magnitude  decay  from  coarse  scales 
to  fine  scales.  The  sampling  patterns  in  the  frequency  transform  domain  are  generated 
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Figure  10:  Optical  sectioning  imaging  result,  (a)  Imaging  target  captured  by  a  CCD  camera, 
which  is  a  microcircuit  chip,  (b)-(g)  Optical  sections  reconstructed  from  our  experiment 
testbed.  The  distance  between  two  adjacent  optical  sections  in  the  depth  direction  is  l//m. 


randomly  according  to  the  underlying  probability  function. 

Equipped  with  fast  transform  algorithms,  the  proposed  sampling  processes  arc  simple, 
fast  and  can  be  easily  implemented.  By  exploiting  the  a  prion  information  of  natural  images, 
the  CS  performance  obtained  with  the  proposed  sampling  method  is  significantly  improved. 
Compared  with  other  sampling  patterns,  such  as  radial  sampling  pattern  and  variable  density 
spiral  that  also  exploit  some  a  pnon  information  of  the  images  [2.  7],  the  proposed  sampling 
patterns  are  not  heuristically  constructed,  but  are  baser!  on  reliable  statistical  models  of 
wavelet  coefficients  and  thus  the  proposed  sampling  patterns  are  analytically  justified. 

4.1  Background  Reconstruction  Algorithm 

4.1.1  Compressed  Sensing  in  a  Transform  Domain 

A  signal  x  £  is  S  sparse  on  some  basis  =  [ip  ,  t/^,  . . . ,  0v]  if  x  can  be  represented 
by  a  linear  combination  of  S  vectors  of  'I7  with  6'  <<C  ,Y.  The  signal  can  be  expressed  as: 
x  =  where  0  is  an  N  x  1  vector  with  only  S  non-zero  entries.  Compressible  signals, 
such  as  natural  images,  can  also  be  well  described  by  such  sparse  signal  model  [9].  Here  we 
consider  compressed  image  sampling  in  the  transform  domain^  .  To  obtain  the  sparse  signal 
information,  we  acquire  a  small  set  of  transform  coefficients  of  x  in  <h.  The  measurements 
are  given  by:  y  =  ‘Fnx,  where  *I>si  is  a  M  x  N  matrix  with  M  <C  N  and  y  =  [j/i ,  y-2,  ■ . . ,  Vm}1 
represents  the  M  measurements.  Each  row  of<I>  </>.,  is  taken  from  a  subset  fi  C  {1, . . . ,  N} 
of  the  atoms  of  <F.  If  and  $  are  incoherent  with  each  other  and  0  is  randomly  chosen,  then 
x  can  be  recovered  from  y  with  high  probability  when  M  satisfies:  M  —  CS  log  N  <§;  N, 
where  C  >  1  is  the  oversampling  factor  [1,  2].  The  incoherence  between  two  bases  (  4$)  is 
characterized  by  the  mutual  coherence  defined  as:  =  sup{|(^,0)|  :i/’€ 'F,  06*1}  [10]. 

It  is  found  that  the  low  coherence  property  holds  for  many  pairs  of  bases  (T,  <I>). 

In  practical  applications  where  noise  is  present,  the  measurements  are  modeled  as:  y„ 
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Figure  11:  Optical  sectioning  imaging  result  of  a  pollen  grain  specimen.  The  distance 
between  adjacent  optical  sections  in  the  depth  direction  is  1  pun. 


Tjj'M  +  n,  where  n  is  zero-mean  additive  white  Gaussian  noise.  The  reconstructed  signal 
can  be  obtained  by  using  the  Basis  Pursuit  Denoising  (BPDN)  algorithm  which  solves  the 
following  problem  [IT: 


6  =  argmin  ||<I>s2'I'0  -  yri|||  +  A||0||, 


0) 


where  ||#||i  =  JT  |#j|  and  A  >  0  depends  on  the  noise  level.  Note  that  BPDN  works 
equally  well  for  the  approximate  reconstruction  of  compressible  signals  [llj.  If  there  exist 
fast  algorithms  associated  with  both  <J>  and'F  ,  then  a  fast  reconstruction  algorithm  can  be 
implemented  for  signal  reconstruction  [4]. 

4.1.2  Incoherent  Fourier  Sampling  of  Sparse  Wavelet  Subbands 

Wavelets  have  well  defined  spectral  characteristics  in  the  Fourier  domain  [121.  A  coarse  scale 
wavelet  has  its  spectrum  localized  in  the  low'  frequency  band  whereas  a  fine  scale  w'avelet 
has  its  spectrum  widely  spread  out  in  the  high  frequency  band.  In  [4],  it.  is  observed  that 
if  the  signal  to  be  acquired  is  restricted  to  a  certain  wavelet  subband,  a  set  of  incoherent 
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measurements  can  be  obtained  by  selecting  the  Fourier  coefficients  from  the  corresponding 
frequency  band  where  significant  energy  of  the  wavelet  exists. 

Let  k  =  1, 2 ,K  denote  the  scale  of  the  ID  wavelets  where  k  =  1  is  the  finest  scale 
and  k  =  K  is  the  coarsest  scale.  Without  loss  of  generality,  it  is  assumed  that  the  wavelet 
has  length  N  =  2h .  Let  ipk,i  denote  the  ID  wavelet  at  a  scale  k  with  a  shift  l  G  [0,  N 2  k  —  1], 
then  the  DFT  spectrum  of  (pkj  is  over  the  band  Bk  =  [— ./V2  k ,  —N2  k  *]  U  [7V2  k  *,  N2  fcj. 
To  reconstruct  <pkjt  from  its  DFT  samples,  the  DFT  atoms  are  randomly  selected  within  the 
band  Bk  [4].  Assume  the  number  of  significant  wavelets  over  Bk  is  Sk,  then  approximately 
2 Sk  to  3 Sk  Fourier  measurements  are  needed.  The  probability  that  a  DFT  atom  is  sampled 
depends  on  the  size  of  Bk  and  on  the  sparsity  Sk-  Since  smooth  signals  have  most  of  their 
significant  wavelet  components  clustered  in  the  coarse  scales,  incoherent  Fourier  subband 
sampling  implies  that  the  low  frequency  Fourier  coefficients,  which  contain  much  of  the 
signal’s  energy,  should  be  sampled  with  higher  probability  1 . 

The  principle  of  Fourier  sampling  of  wavelet  subbands  can  be  readily  applied  to  2D 
images.  Assume  the  underlying  image  is  of  size  N  x  N.  Let  <fk,ij  denote  the  (z,  j)th  2D 
wavelet  at  the  fcth  scale  of  subband  B  where  B  G  LH.  HL,  H  H,  i,j  —  0, . . . ,  Nk  —  1  and 
Nk  =  N/2k.  The  DFT  spectrum  of  all  the  wavelets  at  the  fctli  scale  can  be  characterized 
by  Bk  =  ([— AT2_fc,  —N2  k  *]  U  \N2~k~\N2  *])  x  ([-7V2  k,-N2 U  \N2  k  \N2  '*]). 
Again,  the  probability  that  a  2D  DFT  atom  within  Bk  is  sampled  depends  on  the  size  and 
on  the  number  of  significant  wavelets  of  Bk- 

4.2  Variable  Density  Sampling  in  the  Fourier  Domain 

4.2.1  Variable  Density  Sampling  Functions 

For  2D  natural  images,  the  significant  wavelet  components  tend  to  cluster  around  coarse 
scales.  Thus,  according  to  the  principle  of  Fourier  sampling  of  wavelet  subbands,  the  Fourier 
coefficients  to  be  measured  should  also  cluster  over  low  frequency  bands.  Obviously,  by 
placing  the  samples  selectively  but  still  in  a  random  manner,  we  can  achieve  better  quality 
of  image  reconstruction  than  if  those  measurements  are  chosen  uniformly  random.  In  the 
following,  we  discuss  how  the  sampling  probability  of  the  Fourier  samples  should  change  over 
different  frequency  bands. 

For  natural  images,  it  is  well-known  that  the  distribution  of  wavelet  subband  coefficients 
can  be  adequately  described  by  GGD  models  [13].  Let  wfjk  denote  the  (b  j)th  wavelet 
coefficient  in  the  Arth  scale.  For  ease  of  analysis,  we  assume  that  wf3  k,  for  B  G  LH,  HL .  H  H . 
belongs  to  the  same  GGD  with  the  probability  density  function  given  by: 

f(x]ak,0k)  =  7}(ak,  0)  Qxp{-[\x\/d(ak-  0k)}lh}  -  (2) 

where  d(ok,(3k)  =  ak  ^ r/(<7fc,/3fc)  =  2r(7t.^(gfcA)  ™dr(  *)  =  f(~  t  +  1  dt  is  the 

Gamma  function.  The  parameters  and  /3k  arc  the  standard  deviation  and  the  shape 
parameter  of  the  Arth  scale,  respectively.  Typical  values  of  (3k  for  natural  images  falls  in  the 

Results  in  [4]  do  not  directly  imply  that  if  multiple  subbands  contain  nonzero  elements  and  the  corre¬ 
sponding  frequency  intervals  are  sampled  proportionally  to  their  occupancy,  the  signal  can  st  ill  be  recovered 
with  high  probability.  Thus,  splitting  the  signal  into  various  wavelet  subbands  and  then  performing  inco¬ 
herent  sampling  is  suggested. 
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range  [0.5, 1]  [13].  The  variance  of  the  A:th  scale  wavelet  coefficients  is  given  by  ak  which  is 
assumed  to  decrease  exponentially  from  coarse  scales  to  fine  scales  [14]: 

a\  =  2-q(k-*),  (3) 


where  a  >  0.  It  was  shown  in  [14]  that  a  can  range  from  2.25  to  3.1  based  on  the  inference 
from  empirical  studies.  Equation  (3)  also  reveals  the  fact  that  the  magnitude  of  wavelet 
coefficient  decays  rapidly  from  coarse  scales  to  fine  scales.  Here  the  coarsest  scaling  coefficient 
of  the  image  is  assumed  to  belong  to  a  uniform  distribution  C7 (0, 1). 

We  first  discuss  how  the  sampling  probability  of  DFT  atoms  should  change  according 
to  the  incoherent  sampling  of  wavelet  subbands.  Here  we  define  the  significant  wavelet 
coefficients  as  the  wavelet  coefficients  whose  magnitudes  are  larger  than  a  threshold  /i  >  0. 
Since  the  number  of  wavelet  coefficients  at  scale  k  is  3(4^  k)  ,  the'  mean  of  the  number  of 
significant  wavelet  coefficients  at  scale  k  is: 


A*;  =  6(4a  * 


f[x\ak,Pk) 


dx. 


(4) 


Consider  the  (m,  n)th  DFT  atom  within  the  band  Bk  corresponding  to  the  wavelet  scale 
k  =  K  —  |_log2(2.s)J  where  s  =  max{|m|,|n|}  and  —N/2  <  rn,n  <  N/2.  It  is  clear  that 
the  number  of  DFT  atoms  within  Bk  satisfies:  <;k  oc  4A  k.  Thus,  based  on  the  incoherent 
Fourier  subband  sampling  method,  the  probability  that  the  (m,n)th  DFT  atom  is  selected 

also  satisfies:  ^  r°° 

p(rn,n)  oc  —  oc  /  f(x;<r k,fik)  dx.  (5) 

<A-  ./„ 

Before  examining  p(m,  n)  in  the  general  case  it  is  insightful  to  consider  a  special  case  of 
the  GGD:  the  Laplacian  (f3k  =  1).  The  Laplacian  distribution  is  of  particular  importance  for 
the  following  two  reasons:  1)  it  is  analytically  more  tractable;  2)  the  empirical  distribution 
of  subband  wavelets  for  many  natural  images  are  close  to  Laplacian  [13].  With  f3k  =  1  and 
d(o:fc,/3fc)  =  0*/\/2,  it  can  be  shown  that: 


Thus,  p(m,  n)  can  be  approximated  as: 


\  C  XP 


(G) 


p(m,  n)  oc  cxp[— (2^  pss)]. 


(7) 


It  is  clear  that  p(m,n)  decays  exponentially  alone  with  the  atom  coordinates  s  and  the  decay 
rate  depends  on  the  image  characteristic  parameter  a.  Analysis  based  on  the  Laplacian 
distribution  indicates  that  in  the  general  case,  p(rn,  n)  should  also  decay  exponentially  with 
its  atom  coordinates.  For  the  random  selection  of  DFT  atoms,  it  is  convenient  to  construct 
a  sampling  density  function  in  the  DFT  domain  and  generate  a  sampling  pattern  according 
to  the  sampling  density  function.  Designing  a  sampling  density  function  following  (7)  would 
require  image  dependent  parameter  a,  (3k  and  //.  However,  in  reality,  only  the  number  of 
measurements  J  is  known.  A  method  which  does  not  require  image  dependent  information 
to  design  an  effective  sampling  density  function  is  more  desirable. 

To  conform  to  the  decaying  behavior  of  the  sampling  probability  with  increasing  coordi¬ 
nates  while  keeping  a  simple  form,  we  propose  a  new  family  of  sampling  density  functions 
containing  only  exponential  terms.  Assuming  the  size  of  the  underlying  image  is  M  x  N, 
the  probability  that  the  (m,  n)th  coefficient  is  sampled  is  given  by: 
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PF(m,n)  =  exp 


(8) 


LMSEHSZTl 

” F 

where  —M/2  <  m  <  M/2  and  —N/2  <  n  <  N/2.  aF  >  0  is  a  parameter  characterizing 
the  decay  of  the  sampling  probability.  oF  >  0  is  a  parameter  tuned  to  obtain  the  desired 
number  of  samples.  Since  0  <  pF(m,n)  <  1,  pF(m,n)  is  suitable  as  a  probability  function. 
Furthermore,  the  proposed  sampling  density  function  is  flexible  to  accommodate  various 
decaying  rates  and  different  number  of  measurements  by  timing  the  parameters  aF  and 
aF.  Given  aF,  the  sampling  ratio  defined  as  J/MN  is  a  monotonically  increasing  function 
of  aF.  Thus,  the  search  of  oF  for  a  desired  number  of  measurements  is  straightforward. 
Furthermore,  of  can  also  be  found  by  numerically  solving  the  following  problem: 


(y/x2  +  y2)af 


a 


F 


dxdy  = 


J 


M  x  .V 


(9) 


The  sampling  patterns  generated  from  the  sampling  density  function  are  binary  where  1  at 
(m,n)  indicates  a  sampling  point  and  0  means  no  measurement  on  that  point  is  made.  With 
a  probability  given  by  pF(rn,  u),  1  will  be  generated  at  (ra,  n);  otherwise,  0  will  be  generated. 

Based  on  the  observation  from  the  case  of  the  Laplacian  distribution,  it  is  reasonable 
to  set  aF  =  |a/5,  where  j3  is  defined  as  the  averaged  value  of  several  coarsest  scale  shape 
parameters.  However,  it  should  be  noted  that  the  above  analysis  based  on  the  incoherent 
Fourier  subband  sampling  principle  only  provides  approximately  equal  probability  of  recon¬ 
struction  for  significant  wavelet  components  both  at  coarse  scales  and  at  fines  scales.  For 
image  reconstruction,  the  contribution  of  an  individual  significant  wavelet  component  to 
the  image  quality  should  be  considered.  Equation  (3),  along  with  empirical  data  analysis, 
reveals  that  a  significant  wavelet  coefficient  at  a  coarse  scale  usually  has  much  larger  mag¬ 
nitude  than  that  of  a  significant  wavelet  coefficient  at  a  fine  scale.  Thus,  probabilistically, 
a  coarse  scale  wavelet  has  much  more  contribution  to  the  reconstructed  image  quality  than 
that  of  a  fine  scale  wavelet.  Consequently,  to  achieve  good  reconstructed  image  quality,  it  is 
desirable  to  give  more  preference  to  the  reconstruction  of  coarse  scale  wavelet  components. 
Such  intuitive  analysis  shows  that  an  effective  sampling  density  function  should  decay  faster 
than  predicted  by  aF  —  \a$  and  it  is  desirable  to  increase  the  decaying  parameter  aF  away 
from  Simulation  results  show  that  aF  ~  L3  a  leads  to  good  image  reconstruction,  as 
will  be  illustrated  in  Sec.  4.4.  It  is  also  observed  that  the  quality  of  the  reconstructed  image 
is  robust  to  variations  of  aF  if  aF  is  sufficiently  large.  If  an  estimation  of  a  is  not  available, 
then  setting  aF  —  3.5  is  recommended  since  this  value  leads  to  robust  and  satisfactory  CS 
performance. 


4.2.2  Variable  Density  Sampling  Spiral  for  MRI 

The  variable  density  sampling  pattern  proposed  in  Sec.  4.2.1  can  be  applied  to  MRI  in  k- 
space,  leading  to  a  reduced  set  of  measurements.  However,  this  method  does  not  necessary 
reduce  the  total  scan  time.  For  MRI  applications  that  have  limited  scan  time,  like  functional 
MRI,  a  single-shot  spiral  pattern  can  effectively  cover  the  A>space  within  one  repetition  time 
(Tr)  period  [15]. 

Based  on  the  proposed  variable  density  sampling  function  pF(?n,n),  a  single-shot  spiral 
sampling  pattern  is  designed  which  involves  two  steps.  In  the  first  step,  starting  from  the 
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origin  of  the  fc-space,  the  proposed  spiral  samples  all  the  points  within  radius  7  on  the 
sampling  grid.  For  a  given  threshold  0  <  </,  <  1,  7  is  set  such  that  pF(m,n)  >  tn  for 
\/m2  +  n2  <  7.  I11  the  second  step,  where  \/ rn2  +  n2  >  7,  the  trajectory  of  the  spiral  is 
described  by  the  polar  equation: 

r(i)  =  exp  (&“'),  F  <t<  (-y-)  (10) 

where  b  >  0  and  N  is  the  size  of  the  sampling  grid  in  the  fc-space.  The  parameter  t  is 
sampled  sufficiently  so  that  the  sampling  points  are  continuously  located  along  the  spiral 
and  b  is  adjusted  to  obtain  a  desired  number  of  samples.  It  is  shown  that  th  —  0.8  leads  to 
spirals  whose  performance  approximates  those  patterns  generated  from  pF(m,n).  Here  it  is 
assumed  that  the  MRI  gradient  hardware  is  designed  such  that  the  slew-rate  and  amplitude 
constraints  allow  the  sampling  trajectory  to  move  from  one  sampling  point  to  its  neighboring 
points  on  the  sampling  grid  without  any  difficulty  [16]. 

4.3  Sampling  in  the  ordered  DHT  domain 

Sampling  in  the  ordered  Discrete  Hadamard  Transform  (DHT)  domain  is  suitable  for  image 
sampling  with  hardware  capable  of  representing  binary  measurement  matrices  since  the 
basis  images  of  2D  ordered  DHT  are  binary.  Furthermore,  ordered  DHT  has  fast  transform 
algorithms.  In  ordered  DHT,  the  atoms  are  ordered  by  the  number  of  sign  changes  (zero 
crossing)  between  consecutive  entries  in  a  Hadamard  atom  [8].  Thus,  ordered  DHT  can  be 
regarded  as  a  generalized  class  of  DFT  and  more  specifically,  a  binarized  version  of  DCT 
sharing  thus  many  properties  of  DFT  and  DCT  [17]. 

To  design  the  incoherent  sampling  pattern  in  the  ordered  DHT  domain,  we  need  to 
exploit  the  fact  that  ordered  DHT  can  be  regarded  as  a  binary  approximation  of  the  DCT. 
For  illustrative  purposes,  we  lirst  describe  the  spectrum  characteristics  of  ID  wavelets  in 
the  DCT  domain.  Let  Vkj{m),  0  <  m  <  N  —  1  be  the  DCT  of  a  ID  wavelet  <Pk,i(n)i 
n  =  0, 1, . . . ,  tV  —  1.  Define  the  averaged  DCT  spectrum  of  wavelets  at  scale  k  as: 

JV2  k  1 

V(m)k  =  (11) 

/=0 

It  can  be  shown  that  the  DCT  spectrum  of  wavelets  at  scale  k  has  approximately  the 
same  shape  as  their  Fourier  spectra.  The  DCT  spectrum  band  Bk  of  wavelets  at  scale  k 
is  Bk  ^  [N2~k,  N2~k+1}.  Similar  to  incoherent  subband  sampling  in  the  Fourier  domain,  to 
reconstruct  ipkyi  from  its  DCT  samples,  we  need  to  select  the  DCT  atoms  randomly  within 
the  band  Now  consider  the  incoherent  sampling  in  the  ordered  DHT  domain.  We  can 
define  the  averaged  ordered  DHT  spectrum  of  wavelets  at  scale  k  as: 

N2  '  *  —  1 

H(m)k=  |/ifclj(m)|,  (12) 

(=0 

where  hkj(m ),  0  <  rn  <  N  —  1,  is  the  ordered  DHT  of  <pk,i(n)-  Since  ordered  DHT  is  a  binary 
approximation  of  the  DCT,  H(m)k  can  be  approximated  as:  H(m)k  ~  V(rn)k-  Figure  12 
shows  H{m)k  for  k  =  1  to  4  where  N  =  256  and  the  underlying  wavelets  are  Daubechies-4 
wavelets.  The  corresponding  DCT  spectrum  at  scale  k  is  also  shown.  It  can  be  seen  that 
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Figure  12:  Ordered  DHT  and  DCT  spectra,  of  Daubechies-4  wavelets,  (a)  at  scale  4;  (b)  at 
scale  3;  (e)  at  scale  2;  (d)  at  scale  1. 


the  averaged  ordered  DHT  spectrum  and  the  averaged  DCT  spectrum  are  similar,  which 
indicates  that  the  principles  of  sampling  in  the  DCT  domain  should  be  equally  applied  to 
the  sampling  in  the  ordered  DHT  domain.  The  above  analysis  can  be  readily  applied  to 
2D  cases.  Following  a  similar  procedure  to  that  described  in  See.  4.2.1,  the  variable  density 
sampling  function  in  the  2D  ordered  DHT  domain  is  designed  as  follows: 

{VW+WY" 


PH  ( m » n)  =  <’XP 


'a 


(13) 


where  0  <  m  <  M  —  1,  0  <  n  <  Ar—  1.  oh  depends  on  the  number  of  samples  and  an  depends 
on  the  desired  decay  of  the  sampling  probability.  The  sampling  patterns  are  then  obtained 
as  realizations  of  pH(m,n ).  As  sampling  in  the  DFT  domain,  the  sampling  ratio  is  a  mono- 
tonieally  increasing  function  of  aH  and  searching  of  a H  for  a  desired  number  of  measurements 
is  straightforward.  Similar  to  (9),  an  initial  guess  of  a H  can  also  be  found.  As  in  Fourier 
sampling,  simulation  results  show  that  aH  %  1.3  a  leads  to  good  image  reconstruction.  If  an 
estimation  of  a  is  not  available,  then  setting  aH  =  3.5  is  also  recommended. 


4.4  Simulations  Reconstruction  Algorithm 

In  this  section,  extensive  simulations  are  provided  to  illustrate  the  effectiveness  of  the  pro¬ 
posed  methods.  The  selection  of  the  design  parameters  are  also  discussed.  In  the  simulations, 
all  the  images  are  assumed  sparse  in  the  Daubechies-8  wavelet  domain  and  pixel  values  are 
within  [0,  lj.  It  is  also  assumed  that  each  measurement  is  corrupted  by  additive  white  Gaus¬ 
sian  noise  with  variance  a2  =  le  4.  The  BPDN  algorithm  is  used  for  image  reconstruction. 

In  the  first  example,  the  proposed  variable  density  sampling  is  applied  to  MRI  image 
acquisition  and  reconstruction.  The  MRI  images  to  be  measured,  Brain  and  Angiography , 
are  of  size  256  x  256  and  are  shown  in  Fig.  13(a)  and  Fig.  13(b),  respectively.  The  esti¬ 
mated  parameter  a  for  each  image  is  2.32  and  2.47,  respectively.  The  sampling  patterns 
under  investigation  contains  15000  samples  and  there  are  15000  real-valued  measurements 
made  along  with  each  sampling  pattern.  Figure  13(c)  shows  the  proposed  single-shot  spiral 
sampling  pattern  for  fast  MRI.  The  parameters  of  the  spiral  trajectory  is  set  as:  aF  =  3.5, 
th  =  0.8,  7  =  47  and  b  =  4.5  x  10  11 .  To  generate  the  spiral  pattern  with  continuously 
located  sampling  points,  the  sampling  interval  along  t  is  set  as:  At  =  0.007.  Figure  13(h) 
shows  part  of  the  reconstructed  Brain  image  corresponding  to  the  region  within  the  white 
frame  in  Fig.  13(a).  It  is  clear  that  with  an  undersampling  ratio  of  22.9%,  the  MRI  image 
is  reconstructed  with  only  small  distortion.  Figure  13(d)  shows  an  example  of  the  sam¬ 
pling  pattern  generated  directly  from  the  proposed  sampling  function  given  by  Eq.  (8)  with 
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Figure  13:  The  original  MRI  images:  (a)  Brain ;  (b)  Angiography,  (e)  Proposed  single- 
shot  spiral  sampling  pattern,  (d)  Proposed  variable  density  sampling  pattern,  (e)  Radial 
sampling  pattern,  (f)  Logarithmic  spiral  pattern,  (g)  Random  phase-encodes  undersampling 
pattern.  Each  sampling  pattern  has  a  22.9%  undersampling  ratio.  Part  of  the  reconstructed 
image  with  (h)  proposed  spiral  sampling  pattern;  (i)  proposed  variable  density  sampling 
pattern:  (j)  radial  sampling  pattern;  (k)  logarithmic  spiral  pattern.  (1)  random  phase-encodes 
undersampling  pattern.  See  Table  1  for  corresponding  PSNR. 


aF  =  3.5  and  aF  =  0.134,  whose  performance  serves  as  a  benchmark  for  the  proposed  spiral 
patterns.  The  corresponding  reconstructed  image  is  shown  in  Fig.  13(i),  which  contains  less 
distortion  than  Fig.  13(h). 

To  compare  the  proposed  sampling  patterns  with  other  practical  sampling  patterns,  we 
also  reconstruct  the  MRI  image  from  samples  taken  from  radial  sampling  pattern,  loga¬ 
rithmic  spiral  sampling  pattern  and  random  phase-encodes  undersainpling  pattern  shown  in 
Fig.  13(e),  Fig.  13(f)  and  Fig.  13(g),  respectively.  The  random  phase-encodes  imdersam- 
pling  is  restricted  to  undersampling  of  the  phase-encodes  and  fully  sampled  readouts  and 
the  sampling  density  scales  according  to  a  power  of  distance  [7].  The  corresponding  recon- 
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Table  1:  Reconstruction  of  MRI  images.  The  performance  is  measured  by  PSNR  (dB).  “P. 
VD”  is  the  proposed  variable  density  sampling;  ilP.  SS”  is  the  proposed  single-shot  spiral 
sampling;  “P-E  U.”  is  the  phase-cncodes  undersampling. _ 


No. 

rvlT 

P.  SS 

Radial 

Log-spiral 

r-Trcr 

5000 

21.03 

19.78 

Brain 

16.26 

17.19 

17.15 

10000 

23.86 

22.45 

19.55 

20.58 

20.46 

15000 

25.56 

24.74 

22.46 

23.40 

22.99 

20000 

27.23 

26.51 

24.68 

25.97 

24.73 

25000 

28.77 

28.32 

27.05 

28.10 

27.01 

5000 

23.68 

Anaioqraphy 
22.32  "  19.87 

20.34 

21.27 

10000 

26.02 

25.07 

22.89 

23.85 

23.39 

15000 

27.85 

27.38 

25.30 

26.10 

25.70 

20000 

29.74 

29.34 

27.25 

28.54 

27.73 

25000 

31.34 

31.09 

29.18 

30.39 

29.61 

structed  MRI  images  for  the  above  sampling  patterns  are  shown  in  Fig.  13(j),  Fig.  13(k)  and 
Fig.  13(1),  respectively.  All  reconstructed  MRI  images  show  more  low  frequency  interference 
and  the  image  qualities  are  not  as  good  as  that  of  the  proposed  single-shot  spiral  sampling 
pattern. 

Table  1  summarizes  more  extensive  simulation  results  for  the  MRI  test  images  recon¬ 
struction  with  different  sampling  patterns,  where  the  number  of  measurements  ranges  from 
5000  to  25000  (undersampling  ratio  from  7.63%  to  38.12%).  The  quality  of  the  reconstructed 
images  is  measured  by  the  peak  signal  to  noise  ratio  (PSNR.)  and  the  data  are  collected  from 
an  average  of  ten  runs.  The  results  of  the  sampling  patterns  generated  directly  from  the 
proposed  sampling  functions  serve  as  the  benchmarks  and  among  all  the  sampling  patterns, 
the  best  result  is  highlighted  in  bold  font  in  each  test  set.  The  simulations  again  show  that 
the  proposed  single-shot  spiral  sampling  patterns  consistently  have  better  performance  than 
other  fast  sampling  patterns.  The  performance  gain  is  2  ~  3  dB  at  ./  =  5000  and  0.2  ~  0.7 
dB  at  J  =  25000  for  both  images.  Furthermore,  the  difference  in  performance  between  the 
proposed  single-shot  spiral  patterns  and  the  sampling  patterns  generated  from  the  proposed 
sampling  density  function  is  small. 

In  the  second  example,  the  proposed  sampling  method  is  applied  to  the  acquisition  and 
reconstruction  of  a  set  of  natural  images:  boat ,  Lena,  Goldhill,  baboon  and  Pentagon  with 
size  256  x  256.  The  estimated  parameter  a  for  each  images  is  estimated  as:  2.65,  2.88, 
2.32,  2.64,  2.27  and  2.05,  respectively.  Figure  14(b)  shows  a  sampling  pattern  in  the  ordered 
DHT  domain  that  contains  20000  sampling  points  generated  from  Eqn.  ( 13)  with  aH  =  3.5 
and  aH  =  0.501.  Part  of  the  reconstructed  boat  image  corresponding  to  the  region  within 
the  white  frame  in  Fig.  14(a)  is  shown  in  Fig.  14(e),  which  shows  that  the  reconstruction 
contains  little  distortion. 

For  comparison  purposes,  we  also  reconstruct  the  boat  image  from  the  samples  taken  from 
the  radial  sampling  pattern  and  logarithmic  spiral  sampling  pattern  shown  in  Fig.  14(c)  and 
Fdg.  14(d),  respectively.  The  number  of  measurements  taken  from  both  sampling  patterns  are 
the  same  as  that  of  the  proposed  sampling  pattern.  Part  of  the  corresponding  reconstructed 
images  are  shown  in  Fig.  14(f)  and  Fig.  14(g),  respectively.  Again,  both  reconstructed  images 
show  more  aliasing  artifacts,  which  indicates  that  they  are  not  as  effective  as  the  proposed 
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Figure  14:  (a)  The  original  image  Boat,  (b)  Proposed  variable  density  sampling  pattern,  (c) 
Radial  sampling  pattern,  (d)  Logarithmic  spiral  sampling  pattern.  Each  sampling  pattern 
has  a  30.5%  undersampling  ratio.  Part  of  the  reconstructed  image  with  (e)  proposed  variable 
density  sampling  pattern;  (f)  radial  sampling  pattern;  (g)  logarithmic  spiral  pattern;  (h) 
SBHE.  See  Table  2  for  corresponding  PSNR. 

sampling  pattern. 

Simulation  results  are  summarized  in  Table  2.  where  the  number  of  measurements  ranges 
from  5000  to  25000  and  the  best  results  as  the  average  of  ten  runs  is  also  highlighted.  For 
each  test  image,  it  is  shown  that  the  proposed  variable  density  sampling  achieves  much 
better  performance  than  logarithmic  spiral  patterns  and  radial  patterns.  The  performance 
gain  is  2  ~  4  dB  at  J  =  5000  and  0.5  ~  0.8  dB  at  J  =  25000.  To  verify  that  the  proposed 
sampling  pattern,  which  exploits  the  a  prion  information  of  natural  images,  achieves  better 
performance  than  methods  that  do  not  exploit  the  a  priori  information,  the  simulation  results 
using  the  Noiselet  ensemble,  SFE  and  SBHE  under  different  number  of  measurements  are  also 
presented  in  Table  2.  To  acquire  the  image  information,  samples  of  the  Noiselets,  SFE  and 
SBHE  are  taken  randomly.  The  proposed  sampling  method  achieves  the  best  performance 
in  all  simulations.  For  illustrative  purposes,  Fig.  14(h)  shows  the  reconstructed  result  for 
image  boat  using  SBHE.  Such  comparison  clearly  shows  the  performance  gain  achieved  by 
exploiting  the  a  priori  information. 

Finally,  we  show  that  the  reconstruction  performance  is  not  sensitive  to  the  parameters 
aF  or  aH  in  the  sampling  functions.  Here,  we  test  the  reconstruction  performance  as  the  pa¬ 
rameter  aH  changes.  Figure  15  shows  the  reconstruction  results  of  all  test  images  with  20(X)0 
number  of  measurements.  Data  is  collected  as  the  result  of  20  runs.  Each  test  image  has 
different  curves  and  textures,  thus  has  different  statistical  model  parameters.  The  sampling 


19 


Table  2:  Reconstruction  of  images  in  ordered  DHT  domain.  The  performance  is  measured 
by  PSNR  (dB).  “P.  VD”  is  the  proposed  variable  density  sampling. 


fso. 

PTVD 

Radial 

Log-spiral 

Noiselet 

”sfi r~ 

SBHE 

5000 

10000 

15000 

20000 

25000 

22.75 

24.93 

26.70 

28.36 

29.55 

19.88 

22.93 

25.09 

26.96 

28.79 

boat 

19.79 

22.77 

25.03 

26.95 

28.64 

18.51 

20.78 

22.67 

24.19 

25.58 

19.27 
21.39 
23.25 
24.93 

26.27 

18.62 

20.74 

22.69 

24.56 

26.18 

5000 

10000 

15000 

20000 

25000 

25.51 

27.91 

29.77 

31.24 

32.57 

21.28 

25.33 

28.24 

30.07 

31.82 

Lena 

21.51 

25.49 

28.04 

29.91 

31.54 

20.14 

23.10 

25.06 

26.76 

28.48 

21.25 

23.87 

25.77 

27.56 

29.17 

20.05 

23.10 

25.45 

27.43 

29.01 

5000 

10000 

15000 

20000 

25000 

24.96 

27.09 

28.84 

30.18 

31.26 

22.14 

25.14 
25.82 
27.30 
30.63 

Goldhill 

21.24 

24.80 

26.96 

28.67 

29.88 

20.58 
22.89 

24.59 
25.98 
27.37 

21.31 

23.36 

25.10 

26.59 

27.93 

20.38 

22.81 

24.73 

26.52 

27.94 

5000 

10000 

15000 

20000 

25000 

22.47 

23.32 

24.31 

25.30 

26.40 

21.09 

22.55 

23.46 

24.60 

25.92 

baboon 

20.93 

22.57 

23.64 

24.67 

25.82 

19.59 

21.08 

22.12 

23.04 

23.88 

20.16 

21.36 
22.33 

23.36 
24.31 

19.48 

20.91 

22.01 

23.05 

24.15 

5000 

10000 

15000 

20000 

25000 

24.78 

26.32 

27.61 

28.88 

29.86 

22.59 

24.75 

26.32 

27.79 

28.13 

pentaqon 

22.69 

24.69 
26.19 
27.75 
28.97 

21.37 

22.92 

24.33 

25.32 

26.22 

21.95 

23.36 

24.54 

25.66 

26.87 

21  32 
22.79 
24.17 
25.50 
26.77 

Figure  15:  Reconstruction  results  of  images  from  proposed  sampling  patterns  generated  by 
different  aH  with  20000  number  of  measurements. 


patterns  are  constructed  with  aH  ranging  from  0.5  to  6.  It  shows  that  the  sampling  patterns 
lead  to  similar  performance  for  all  the  images  when  aH  ranges  from  2.5  to  4.5.  The  differ¬ 
ence  of  performance  measured  by  PSNR  is  within  0.4  dB.  Thus,  the  image  reconstruction 
is  robust  to  variations  over  nH .  An  empirical  rule  for  the  selection  of  aH  is  aH  =  1.3  a.  If 
the  estimation  of  a  is  not  available,  then  set  aH  =  3.5  is  recommended  since  aH  —  3.5  lies 


20 


in  the  middle  of  the  flat  region  shown  in  Fig.  15.  Note  that  the  reconstruction  performance 
tends  to  become  worse  as  an  increases.  Since  larger  aH  means  more  low  frequency  samples 
are  taken,  the  simulation  shows  that  sampling  only  low  frequency  components  is  not  a  good 
strategy. 

5  Conclusion 

In  this  effort,  we  demonstrated  the  experimental  realizations  of  different  CS- based  imaging 
systems.  We  first  built  a  single  pixel  camera  which  captures  2-dimensional  (2D)  optical 
images  using  a  single  pixel  detector.  Based  on  the  single  pixel  camera  architecture,  we 
built  a  CS-MSI  system,  which  uses  three  beamsplitters  and  four  single  pixel  detectors  to 
acquire  optical  images  at  four  spectral  channels  concurrently,  including  a  red  (650  iim), 
green  (550  nm),  blue  (450  nra),  and  an  NIR  spectral  channel  We  also  built  a  CS-based 
optical  sectioning  microscope  (CSM),  which  used  a  single  pixel  detector  to  realize  the  optical 
sectioning  imaging  performance  without  any  mechanical  scanning.  The  optical  sectioning 
imaging  performance  of  the  CSM  system  was  demonstrated  using  a  reflective  imaging  target 
and  also  a  biological  fluorescent  specimen.  A  family  of  variable  density  sampling  patterns 
are  proposed  for  compressed  sensing  of  natural  images  in  the  DFT,  DCT  and  ordered  DHT 
domain,  which  are  based  on  a  reliable  statistical  model  of  natural  images  in  the  sparse 
wavelet  domain.  The  proposed  sampling  method  is  simple,  fast  and  suitable  for  a  wide 
range  of  applications  such  as  in  [18,  19].  Furthermore,  the  a  priori  information  needed 
is  general  and  no  extensive  data  training  is  needed.  Simulations  show  that  the  proposed 
sampling  scheme  is  effective  for  compressed  sensing  of  images  conforming  to  the  proposed 
wavelet  model. 
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